# sRNA-Seq workflow
## parameters need to be changed

threads=$1
sample_list=$2

last_index=`wc -l ${sample_list}|awk '{print $1}' -`
number_of_files=$[$last_index-1]
refseq=/cluster/group/langzhaobo/zblang/genome/SL4.0/S_lycopersicum_chromosomes.4.00.fa
strna_index=/cluster/group/langzhaobo/zblang/genome/SL4.0/ncRNA_bowtie_index/ncRNA_bowtie_index
genome_index=/cluster/group/langzhaobo/zblang/genome/SL4.0/SL4.0_bowtie_index/S_lycopersicum_chromosomes.4.00
universal_adapter="AGATCGGAAGAG"
small_RNA_5_adapter="GATCGTCGGACT"

for id in `seq 1 $last_index`
do
old_sample=`awk 'BEGIN{FS=","}NR=="'$id'"{print $1}' $sample_list`
sample=`awk 'BEGIN{FS=","}NR=="'$id'"{print $2}' $sample_list`
path=`awk 'BEGIN{FS=","}NR=="'$id'"{print $3}' $sample_list`
newpath=`awk 'BEGIN{FS=","}NR=="'$id'"{print $4}' $sample_list`
head="#!/bin/bash\\n#SBATCH -J ${sample}\\n#SBATCH -p batch\\n#SBATCH -N 1\\n#SBATCH --output=${sample}_log.%j.out\\n#SBATCH --error=${sample}_log.%j.err\\n#SBATCH --cpus-per-task=$threads\\n"

## 1. rename

if [ ! -f $newpath ]
then
mkdir -p $newpath
fi

# input:
gz_1=`ls $path|grep _R1.fastq.gz|sed "s#^#$path\/#"`
gz_2=`ls $path|grep _R2.fastq.gz|sed "s#^#$path\/#"`

# output:
fq_1=$newpath/${sample}_R1.fastq.gz
fq_2=$newpath/${sample}_R2.fastq.gz

# shell
cat1="cat $gz_1 > $fq_1\\n"
cat2="cat $gz_2 > $fq_2\\n"
# 2.fastqc.raw

#shell:
fastqc_raw="fastqc --quiet -o $newpath $fq_1 $fq_2\\n"

# 3.trim

#output
R1_trimmed=$newpath/${sample}_R1.trimmed.fastq.gz
R2_trimmed=$newpath/${sample}_R2.trimmed.fastq.gz
# params:
minlength=14
maxlength=35
#shell:
trim="cutadapt -a $universal_adapter -A $small_RNA_5_adapter -o $R1_trimmed -p $R2_trimmed $fq_1 $fq_2 --minimum-length $minlength --maximum-length $maxlength --discard-untrimmed \\n"


## 3.5 gunzip
R1_unzip="gunzip $R1_trimmed\\n"
R2_unzip="gunzip $R2_trimmed\\n"
## 4.fastqc clean
#shell:
fastqc_clean="fastqc --quiet -o $path $R1_trimmed $R2_trimmed\\n"

## 5. mapping
bam=$newpath/${sample}.bam
mapping="bowtie -p $threads $genome_index ${R1_trimmed%%.gz} -k 10 -m 10 -v 0 -S|samtools view -F 4 -bh -  > $bam\\n"

## 6. mapping to structural RNA
stRNA_bam=$newpath/${sample}.stRNA.bam
stRNA_mapping="bowtie -p $threads $strna_index ${R1_trimmed%%.gz} -k 10 -m 10 -v 0 -S|samtools view -F 4 -bh -  > $stRNA_bam\\n"

## 7. extract stRNA reads 
stRNA_reads=$newpath/${sample}.stReadsID.txt
extract_stRNA_reads="samtools view $stRNA_bam |cut -f 1|sort|uniq > $stRNA_reads\\n"

## 8. remove stRNA
rmSt_bam=$newpath/${sample}.rmStruc.bam
rmstruc_cmd="filterbyname.sh in=$bam out=$rmSt_bam names=$stRNA_reads include=false\\n"

## 9. sort bam
rmSt_sort_bam=$newpath/${sample}.rmStruc.sort.bam
sort_bam_cmd="samtools sort $rmSt_bam > $rmSt_sort_bam\\n"


## 10. bam len

#echo -e $head$cat1$cat2$fastqc_raw$trim$fastqc_clean$mapping$stRNA_mapping$extract_stRNA_reads$rmstruc_cmd$sort_bam_cmd > ${sample}.sRNA.pbs
echo -e $head$mapping$stRNA_mapping$extract_stRNA_reads$rmstruc_cmd$sort_bam_cmd > ${sample}.sRNA.pbs
#qsub ${sample}.sRNA.pbs



done